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1 Abstract 

A semi-empirical formalism based on the second moment tight binding approach, consider- 
ing two bands is presented for deriving interatomic potentials for magnetic d-band materials. 
It incorporates an empirical local exchange interaction, which accounts for magnetic effects 
without increasing the computing time required for force evaluation. The consequences of 
applying a two-band picture to transition metal alloys and transition metal impurities is ex- 
amined, which combined with evidence from ab initio calculations leads to some surprisingly 
simplifying conclusions. 

PACS 61.80. Az ll.20.Bc 61.82.Bg 71.15.-m 

2 Introduction 

Semiempirical models for metallic binding have had a long and successful history in computer 
modelling. The distinguishing features of models for metallic bonding are that they are short- 
ranged and non-pairwise. The first of these features arises from the strong screening of the 
nuclear charge by the mobile electrons, the second arises from the delocalisation of those 
electrons. 

The most significant development in accounting for many body effect came in the mid 
eighties with the implementation of 'embedded atom' potentials (EAM)@] (based loosely 
on density functional theory [5]) and 'Finnis-Sinclair' potentials (FS) 6 (based on the tight 
binding second moment approximation 7 ) . The two models have very similar computational 
requirements, and the names are often used interchangably, however there are some distinc- 
tions which come to the fore when considering multicomponent alloys. 

To highlight the differences, the energy according to the EAM is written: 

ueam = Yl v " fa fa ( R v )] w 

ij i j 

Where i and j label atoms of element I and J respectively, V is a pairwise potential 
which depends on both species, Fi and <f>j are the embedding function and charge density 
which depend on one species only. 

In slight contrast, the FS approach implies: 

u FS = £ Vuinj) + £ f[J2 MRij)} (2) 

ij i j 

where <pu is the squared hopping integral between the atoms on site i and site j and F is 
independent of the atomic species {F(x) — y/x in the second moment tight binding approxi- 
mation) In relating theses potential to tight-binding, there is a subtle issue of interpretation 



whether the many body part is a bond or band energyJS] and extensions beyond second mo- 
ment work explicitly with bonds|I|. While mindful of the distinction, for convenience the 
discussion here is presented in terms of bands projected onto atoms. 

This subtle distinction becomes relevant for alloys: under EAM one requires no further 
refitting of the many body function for each species, while the FS implies that a new (f>jj 
function should be fitted for each pair of elements. 

The success of ab initio methods in describing materials behaviour, and the possibility 
of using it directly to parameterise kinetic Monte Carlo simulations, may appear to negate 
the need for empirical potentials. However the lengthscales of correlated events involved in 
certain types of calculation - such as radiation damage cascades from fusion neutrons - are 
so large that empirical potentials still have a role to play. 

ab initio calculation has revealed a number of shortcomings in previous parameterisations 
of interatomic potentials for radiation damage application. In particular, when the original 
EAM parameterisation were done in the 1980s, the only information which affected the fit 
of the potential at short range came from isotropic compression: extrapolation of fits to 
empirical equations of state such as Murnaghan or Rose Jl , or simple quantum electron gas 
calculations % . Recent total energy calculations reveal a picture in which isotropic compres- 
sion of metals is resisted by many-body effects (such as electron kinetic energy) , and that ions 
can approach one another rather more closely than had been predicted. A consequence of this 
is typically that interstitial formation energies are lower than in previous parameterisations 
had predicted, and the associated strain fileds are smaller [3] . 

A philosophical issue arises at this stage. It is clear that such many body potentials 
cannot contain all the correct physics of the systems they are describing. What then should 
they aim to achieve? 

In early work, it was attractive to show how the physical anomalies of pair potentials 
(zero Cauchy pressure, equivalence of vacancy formation and cohesive energy) can be elimi- 
nated with relatively little computational cost. In nuclear materials simulation however, the 
emphasis is less on physical elegance and more on practical application. In the multiscale 
modelling framework molecular dynamics modelling is simply an interpolation tool between 
exact ab initio and experimental results and the non-equilibrium behaviour of many-particle 
systems. Thus there is little justification for constraining the ability to precisely reproduce 
defect properties on the grounds of physical elegance. Consequently, recent parameterisations 
of "many-body" (MB) models write the energy in the most general form: 



which for an N-component systems gives N(N + 1) functions available for fitting. 

Most of the physical discussion motivating the potential centres on the form of the many 
body term, the pairwise part representing "everything else" . However, for high energy colli- 
sions and interstitials it is the pairwise part which tends to control behaviour. 

It is possible to generate still more complex potentials with a similar computational cost 
to EAM. Examples include the MEAM 10 , which incorporates angular effects in a way which 
uses only quantities already calculated in central-force molecular dynamics models, and the 
two-band model which describes the electronic state of each atom in an analytically solvable 
formfH]. 

One important system which is not easily described by these methods is iron. Although 
there have been a number of parameterizations for iron in a single phase, the essential physics 
leading to the stabilisation of the bcc phase is missing. Ab initio calculations show that 
without ferromagnetism, the fee structure is favoured, and that the transition to the high 
pressure (hep) or high temperature (fee) phases is accompanied by a loss of ferromagnetism. 
Similarly, the atoms close to interstitial defects have a reduced magnetisation QJ, which has 
the effect of making them more compressible. 
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In a recent paper [14] . a simple two-band picture was derived for Cs to capture s-d elec- 
tron transfer. This described the second-moment tight binding energy in terms of a single 
parameter (the transfer of electrons between bands) at each site. Surprisingly, the parameter- 
ization proved transferrable across the transition metal series. Crucially for use in molecular 
dynamics, the minimum-energy value of this parameter can be evaluated at each atom in- 
dependently and analytically so the potential has the same computer requirement as EAM. 
Moreover, the energy is variational in this parameter, meaning that forces take a very simple 
form. 

Here it is shown how a similar two-band picture within the second-moment approximation 
can be developed to incorporate band ferromagnetism. In the simplest approximation, it 
give the surprising result that two-band effect can be captured in a standard many-body 
parameterisation, only the interpretation of the embedding function being changed. 

3 Two d-band model - magnetism 
3.1 Modified Many-body term 

In this section we show that incorporating band magnetism in the FS picture while ignor- 
ing the repulsion effect of enhanced Pauli exclusion gives rise to an embedded atom-type 
formalism with a non-monotonic embedding function. 

In the second moment approximation to tight binding, the cohesive energy is proportional 
to the square root of the bandwidth, which can be approximated as a sum of pairwise 
potentials representing squared hopping integrals 6 . Assuming atomic charge neutrality, this 
argument can be extended to all band occupancies and shapes JSj- For simplicity, consider 
a rectangular d-band of full width W centred on Eq . The bond energy for a single spin band 
relative to the free atom is given by: 

r E f =(Z/N-^)w 7 ] 7 \ 

tf= E/WdE= — ( —-l)W/2 (4) 

J-W/2 iV iV 

where Z is the occupation of the band and uparrows denote "spin up" . 

To describe the ferromagnetic case, it is assumed that there are two independent (i-bands 
corresponding to opposite spins, and that these can be projected onto an atom. In the atomic 
case Hund's rules lead to a spin of S=2 for iron, and there is an energy U x associated with 
transferring an electron to a lower spin state. In the solid, the simplest method is to set 
U x to be proportional to the spin with the coefficient of proportionality being an adjustable 
parameter, E . 

U x = -E \Zl - Z±\ (5) 

Defining the spin, S = — Z- 1 and noticing that with charge neutrality the total number 
of electrons at a site is conserved T — (Z^ + Z- 1 ), we make the assumption that the bandwidth 
W is the same for each band (see figC]). The two-band binding energy [Tl| for a d-band with 
capacity N=5 on a site i is then: 

Ui = Uj + Ui + Uf = +^(T 2 + S?) - TWi/2 - E Q S t (6) 

This equation gives us the optimal value for the magnetisation of a given atom of Si = 
2NEo/Wi, and the many-body energy of an atom with T=6, N=5 (suppressing the i label) 
as: 



U= -6W/5-5E$/W E /Wi< 0.4 

-2W/5 - AEq E /Wi > 0.4 (7) 
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Figure 1: BCC iron density of states for majority and minority spin bands from ab initio 
spin-dependent GGA pseudopotential calculations with 4913 k=points, adjusted so the ths 
Fermi energy lies at the zero of energy. The two band model assumes that the bands have 
the same shape and width, but are displaced in energy relative to one another. Here it can 
be seen that this is the case. 

Which introduces the constraint that there is a maximum value for S. For a material 
with T d-electrons (where T > 5) . transfer of electrons between the spin bands becomes 
advantageous for W > 10-E'o/(10 — T). For smaller W the spin | band is full and the energy 
is simply proportional to the bandwidth of the j, band as in the Finnis- Sinclair model. 

Within the second moment model, the bandwidth W is given by the square root of the 
sum of the squares of the hopping integrals[7|. This sum can be represented by a pairwise 
potential and is the same for each band. 



Interatomic potentials also include the effect of the ions and non- valence electrons: this is 
usually represented by a simple pairwise potential. This gives us the final functional form for 
incorporating magnetism in the ferromagnetic regime within the Finnis-Sinclair formalism: 



where H is the Heaviside step function, B is a constant, pj — 4>(fij) , ^ an( i & are 
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empirically fitted pairwise potentials, and the zero of energy is shifted to correspond to the 
non-magnetic atom. 

Note that this form does not explicitly include S, and that it has the form of the embedded 
atom model with an embedding function Fi(x) = \/x(l — B/x). In particular, the local 
magnetisation means that the embedding function now depends on the type of atom through 
the amount of the particular band projected onto it. 

Although this model incorporates magnetism, and provides a way to calculate the mag- 
netic moment at each site, it is possible to use it without actually calculating S - the imple- 
mentation with an additional many-body repulsive term is similar to the many-body potential 
method of Mendelev et al llj. In terms of forces and energetics, it will perform in a similar 
way. In particular, there is no first order magnetic transition. 

The anomaly with this model is that the linear expression for U leads to a non-zero 
magnetisation at all densities (although this becomes vanishingly small at high densities). If 
instead of equationGOwe use a Stoner-type quadratic term U x = IoS 2 this changes the eqnEl 
such that the spin state flips from fully magnetised to fully demagnetised at Wi = 4NIq. 

Neither of these behaviours is quite correct, and a better treatment involves making the 
pairwise term spin-dependent, which is the subject of the following section. 

3.2 Incorporating Pauli Repulsion 

Pauli repulsion arises from electron eigenstates being orthogonal. While its nature on a single 
atom is complex, its interatomic effects can be modelled as a pairwise effect of repulsion 
between electrons of similar spins. The secondary effect of magnetisation is that there are 
more electrons in one band than another, more same-spin electron pairs to repel one another, 
and so the repulsion between those bands is enhanced. Using the Stonor-type U x = IoS 2 , 
the many body energy on an atom becomes 

U = U ] ~ + U l + U x = +j^(T 2 + S 2 ) - TW/2 + I S 2 (10) 

This shows that a high magnetisation changes the energy in two ways, firstly by effectively 
widening the bands and secondly via Io- 

We now consider the pairwise repulsive part of the potential. Conceptually, this contains 
two effects, the standard paramagnetic repulsion which represents the screened Coulomb 
repulsion of the ions and the core-core repulsion, and an additional Si dependent term arising 
from Pauli repulsion between like-spin electrons 1 . Writing this in a separable form: 

V(r tJ ) = V ( rij ) + (Si + S 3 )V m (r l3 ) (11) 
The total energy on the i atom is then 



-TWi/2 + (S 2 + T 2 )Wi/AN + SJ + \ £ V (r tJ ) + S t ^ V m {r l3 ) 



(12) 



The key to this step how we "assign" the Pauli energy to each atom. If we choose to do it 
such that the energy assigned to each atom i depends only on the value of Si at that atom, 
this allows us to minimise Utot with respect to {Si} at each site independently which gives 
us: 



1 Note there are issues about antiferromagnetism here. An AFM state with Si + Sj =0 would have a lower 
repulsive energy. In the tight-binding picture this would be compensated by a much reduced hopping integral 
and hence lower W. If we insist on Si > then we suppress these solutions and can model ferromagnetic or 
diamagnetic iron. Also, as with DFT-GGS/LDA the spin is Ising-like. 
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Values of Si are bounded by and 2N-T. Thus, we can solve analytically for S, and the 
computational cost of the calculation becomes equivalent to a standard EAM. 

The energy depends variationally on the values of 5, : dU / dSi = 0. This allows application 
of the Hellmann-Feynman theorum in the force calculation to eliminate terms involving 
changes in Si. 

Thus the complexity of the forces is also just like EAM 



1 dVo dV m , 



2 I rp2 



4iV 2 / ^ dr„ 



Sj + T 2 T 



AN 2 1 dnj 



(16) 



To make a potential, one needs to fit the pairwise functions V(ry) and (^(fij) an d the 
parameter Iq . This model can be generalised from the Finnis-Sinclair form to the EAM form 
by using a parameterised embedding function rather than the square root. 



4 Two band model for alloys 

4.1 Transition metal alloys 

An alternative use of a two d-band model occurs when an alloy of two d-metals is formed. 
The difficulty with one band on each atom is that charge transfer between sites becomes 
non-local and minimising energy becomes non-analytic J4| . It is possible to use a fictitious 
dynamics for the charge transfer, in the same spirit as Car-Parinello, but this would introduce 
a large number of additional equations of motion. 

The convenient solution, implicit in equation |3 is to assume that a similar d— band can 
be projected onto both sites, shifted for charge neutrality. This would allow us to use the 
MB approach. Here we investigate whether such an assumption can be justified. Fig|5]shows 
a density of states for a MoV alloy. The notable feature is that the alloy bands are narrower 
than the pure metal bands, suggesting that 24>AB{ r ij) < <pAA {r%j ) + ^ssC^ij) 

Furthermore, calculations of the projected density of individual electron states onto in- 
dividual atoms for isolated Mo and V and for Cu in Fe shows that there are no localised 
G?-states - the band structure is projected fairly evenly onto each species. 

4.2 p-band impurities - C, P, Al 

The presence of non-transition metal impurities in steels is of major interest for a wide range 
of applications. The combination of p- and d-based materials makes for a complicated band 
structure and the picture of simply adding electrons to a rigid band is no longer tenable. 
Fig0 shows the ab initio calculated density of states for FesC cementite, compared with 
the density of states calculated for the identical arrangement of Fe atoms with the carbons 
removed. 

The first striking feature is that the occupied carbon bands can be readily distinguished in 
FeaC. Unoccupied carbon p-bands lie well above the Fermi energy, meeting the requirement 
for charge neutrality which underlies the idea of and band energy determined from a local 
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Figure 2: Density of states from ab initio GGA pseudopotential calculations with 4913 
k=points, adjusted so that the Fermi energy lies at the zero of energy. Upper bands: pure 
V (dots) and Mo (solid) at their equilibrium structure (bcc, 3.006Aand 3.184Arespectively). 
Lower bands (plotted upside down for ease of comparison): MoV alloy in CsCl structure 
(lattice parameter 3.086A), with density of states projected onto the two atoms according 
to the integrated electron density within 1.5Aof the nucleus. Further 16-atom supercell cal- 
culations for single substitutional Mo in V and Cu in Fe show similarly undifferentiated 
bands, although when semi-core p-electrons are included in the valance band this method 
show near-perfect localisation. 

density of states [T5]. This suggests that, in the FS picture at least, the hopping integral 
between bands can be ignored. Furthermore, the carbon bands lie far below the Fermi 
energy. This suggests that they will not contribute to the metallic bonding and cannot be 
treated in the same way as the ci-electrons. This means that to a first approximation one can 
set 4>Fec(f) = 0, and treat the iron-carbon interaction purely as a pair potential. 

Similar unhybridised carbon states have been seen for isolated defects: the one occasion 
when the bonding of the dissolved carbon changes qualitatively being the formation of a 
covalent bond in a C dimer located in a Fe vacancy 13 . 

A similar same situation applies to embedding a phosphorus atom in iron. The phosphorus 
s-electrons play no part in the bonding, and the p-states lie below the (i-band (Fig^J. 

5 Multiple elements without refitting 

Finnis-Sinclair type models have had reasonable success in modelling interatomic alloys of 
transition metals, here the extension of the theory is straightforward. If accurate results 
for radiation damage are required, then some refitting is necessary. For alloys with low 
concentrations of impurities, such as reactor steels, a more straightforward approach may 
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Figure 3: Density of states from ab initio GGA pseudopotential calculation with 1000 
k=points. Upper bands: FeaC cementite with atomic positions and unit cell are fully re- 
laxed (a=5.069A, b=6.744A, c= 4.514.A). Lower bands: unit cell containing Fe in the same 
locations as cementite, but no carbon, this gives a non-ferromagnetic ground state. 



be appropriate. With the simple two-band iron model above, we can assume that the main 
effect of small concentrations of local impurities is not to change the shape of the projected 
band structure, but to change the number of electrons in the band. As before, the MB energy 
can be separated into contributions from each atom: 

i 

following the procedure adopted for Cs|14| an approximate potential for an embedded 
transition metal impurity can be made using the same parameters as for iron with lengths 
scaled by the Fermi vector (T -1 / 3 ) and energies by T 1 / 2 . The separation of the pair potential 
into the form: 



V IJ (r ij ) = V I (r ij ) + Vj(r ij ) (18) 

implies that the screening for the impurity is due to the iron charge density. Hence such 
a pairwise potential would not be transferrable to environments of high concentration of the 
impurity atom. It is already well known that potentials fitted to give a good description of 
the dilute environment show poor transferrability across the phase diagram. [TCI IT7] 



8 




Figure 4: Density of states from ab initio GGA pseudopotential calculations with 729 
k=points, adjusted so that the Fermi energy lies at the zero of energy. Upper bands: su- 
percell of 15 iron atoms with a substitutional P impurity, relaxed to its equilibrium volume 
(186A 3 ) with ions relaxed and magnetic moment 35.3/is. Lower bands :supcrcell of 15 iron 
atoms with a vacancy, relaxed to its equilibrium volume (183. 7A 3 ) with ions relaxed and 
magnetic moment 36A^b- Note that the two densities of state are extremely similar, the 
main effect of the P being the addition of two low lying s-states and some additional p-states 
and the low energy end of the <i-band. 



6 Potentials for molecular dynamics - or vice versa 

The concept of designing interatomic potentials is typically to start with a full quantum 
mechanical treatment of the electrons, and then via a series of approximations arive at a form 
suitable for use in molecular dynamics. Hence EAM is formulated as a extreme localised form 
of the density functional theory, while Finnis Sinclair is a simplified tight binding. It is unclear 
to what extent the correctness of the original physics is retained - certainly in applications 
empirical refitting is used to "fold back" all the missing physics in an approximate way. 

The ultimate test of the usefulness of a potential is its transferrability - does it work in 
environments different from where it was fitted? If it is transferrable, it is quite acceptible 
for a potential to include arbitrarily introduced terms to compensate for the uncontrolled 
approximations made from the full quantum mechanical treatment. 

Hence an alternate view of potentials is to start with the needs of molecular dynamics 
and ask what functional forms are possible. The constraints arc that we require short-ranged 
interactions based on atomic positions only. The potential may then use any information 
available from the molecular dynamics to reproduce the forces correctly. In fact, much 
information is available in a real molecular dynamics calculation which might be useful in 
formulating the potential. 
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For simplicity, we assume the calculation does not explicitly evaluate three-body terms 
- this involves an additional calculational overhead and there are no such terms in the full 
many-body quantum mechanical Hamiltonian. The following quantities can be evaluated: 

a) functions of separations at each atom (EAM's "electron density" or FS's "sum of 
squared hopping integrals"). 

b) functions of other rotationally invariant quantities at each atom (as in MEAM) 10 . 

c) on-site parameters with respect to which the energy can be minimised locally and 
analytically (as in the two-band model[Tl]'). 

d) the mean electron density, and the free electron Fermi energy. 



6.1 Evaluating and Using the Fermi energy 

Almost all quantum treatments of metals require calculation of the Fermi energy. In the 
Finnis Sinclair approach it is the fixing of the Fermi level combined with the condition 
of charge neutrality which allows the method to work for other than half-filled bands. [T5] 
Further, it has been known since Hume-Rothery that Fermi surface effects lies at the heart 
of the structural properties of metallic alloys. 

In principle, a short-ranged potential cannot incorporate such effects, because the Fermi 
energy depends on delocalised electrons: an extremely long range effect. However, in setting 
up a molecular dynamics simulation, all the information needed to calculate the Fermi level 
is often included. Most radiation damage calculations entail a constant volume ensemble 
simulation of a bulk material with a few defects: the mean electron density is known and 
constant throughout, and the Fermi level can therefore be established. In surface calculations 
there appears to be a problem, however even here the molecular dynamics uses boundary 
conditions (periodic slab geometry, fixed layer, etc. ) which are intended to represent contact 
with a bulk of fixed density. Again, the Fermi Energy for this can readily be established. In 
constant pressure simulations, it may appear that the Fermi energy should vary, yet even this 
is unclear - in standard Parrinello-Rahman dynamics an expansion of the region of interest 
is accompanied by and expansion of all periodic-image supercells. In the case of a phase 
transition, this is reasonable and the Fermi energy can track the density, however in other 
cases the local expansion would actually cause a compression in the surrounding region, 
preserving the Fermi energy at its bulk value. Hence either treatment of the Fermi level has 
equally good physical justification as the Parrinello-Rahman scheme itself. 

Thus there should be no objection to using the Fermi Energy as part of the formulation of 
a potential, since its value can readily be determined at the outset of a molecular dynamics 
calculation. 

One anomaly this removes is the series of phase transitions under pressure exhibited by 
empirical potentials: e.g. even the OK ground state of the Lennard- Jones potential switches 
between fee and hep as a function of pressure[2I], though no real material does this. Other 
close-packed metal potentials behave similarly. The reason for this is that phase stability 
depends on the long-range part of the potential - beyond second neighbours at equilibrium. 
Under compression various shells of neighbours come into the range of the potential, increas- 
ing or decreasing the energy of the two phases. 

The resolution to this problem has long been known. Thirty-five years ago Heine and 
WeairefTS} showed that within pseudopotential theory phase stability for simple metals the 
energy could be calculated as a volume term independent of atomic positions, plus a pair 
potential with asymptotic form: 



V{R) oc -):7 ^,' (19) 



cos(2kpR 
(2k F R) 3 

The key point here is that kpR is close to constant for isotropic compression. Hence this 
term will not contribute to the change in relative phase stability under pressure. The r -3 form 
can be made usable in molecular dynamics if the long range pairwise part is damped with an 
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cxponential|19|: this corresponds to a finite electron temperature and does not significantly 
affect the forces. The same fictitious-electron-temperature approximation is used routinely 
in ab initio calculations where it is known as " Fermi smearing" . 

Recently, potentials which attempt to measure the density pi with local sampling have 
appeared[2Dl : the most effective local sampling measure seems to be a sum of a Gaussian 
potential. These allow the pair potential to depend explicitly on the electron density and 
take the form: 



However, as discussed above, for practical molecular dynamics such short range approxi- 
mations to find the Fermi energy (and their associated calculational cost) are unnecessary. 

7 Conclusions 

In this paper we have laid out how the ideas underlying the two-band model, already applied 
to caesium|14|. can be extended to magnetic elements such as iron and its alloys. By way of 
ab initio simulation we have shown how the site-projected density of states, which forms the 
conceptual framework behind Finnis-Sinclair type potentials, suggest particular approaches 
to different systems. In particular, if magnetism is introduced without affecting the pairwise 
potential, an EAM-type model is recovered - the onsite magnetisation need only be calculated 
if different Pauli repulsion for majority and minority bands is assumed. Similarly, for sp- 
bonded impurities C and P in iron, we find valence states well below the Fermi energy and 
little effect on the d-bands. This suggests that these impurities can be sensibly treated with 
pair potentials + unmodified many-body iron potential, unless they get close enough to bond 
to one another. 

Accurate parameterisations are not presented here, rather the paper points the way to de- 
veloping future interatomic potentials and provides some post facto justification for previous 
studies. 
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